function relax(i,j)

    global u v w du dv dw a b c dt h m

    
    %
    %  u runs ->
    %  v runs ^
    %
    
    uij = u(i,j);
    vij = v(i,j);
    wij = w(i,j);
    
    %phi_right = sqrt((h + u(i,j+1) - uij)^2 + (vij - v(i,j+1))^2);
    %phi_left = sqrt((h + uij - u(i,j-1))^2 + (vij - v(i,j-1))^2);
    %phi_up = sqrt((h + v(i-1,j) - vij)^2 + (uij - u(i-1,j))^2);
    %phi_down = sqrt((h + vij - v(i+1,j))^2 + (uij - u(i+1,j))^2);
        
    ddu = a*uij + b*du(i,j) ...
        - c*(h - u(i,j-1) + uij) ...    % fu_left
        + c*(h + u(i,j+1) - uij) ...    % fu_right
        + c*(u(i-1,j) - uij) ...        % fu_up
        + c*(u(i+1,j) - uij);           % fu_down

    du(i,j) = du(i,j) + ddu*dt;
    
    ddv = a*vij + b*dv(i,j) ...
        + c*(v(i,j-1) - vij) ...        % fv_left
        + c*(v(i,j+1) - vij) ...        % fv_right
        + c*(h + v(i-1,j) - vij) ...    % fv_up
        - c*(h - v(i+1,j) + vij);       % fv_down
    
    dv(i,j) = dv(i,j) + ddv*dt;    

    ddw = a*wij + b*dw(i,j) ...
        + c*(w(i,j-1) - wij) ...    % fu_left
        + c*(w(i,j+1) - wij) ...    % fu_right
        + c*(w(i-1,j) - wij) ...        % fu_up
        + c*(w(i+1,j) - wij);           % fu_down

    dw(i,j) = dw(i,j) + ddw*dt;    

   % ep(i,j) = 1/2*c*uij^2 + 1/2*c*vij^2;
   % ek(i,j) = 1/2*m*du(i,j)^2 + 1/2*m*dv(i,j)^2;
end